GetRMSEfloat Function

FUNCTION GetRMSEfloat &
!
(grid1, grid2, maskReal, maskInteger, nrmse) &
!
RESULT (rmse)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1
TYPE (grid_real), INTENT(IN) :: grid2
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger
LOGICAL, OPTIONAL, INTENT(IN) :: nrmse


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: rmse, mean
INTEGER             :: countCells
!---------------------------end of declarations--------------------------------

rmse = 0.
countCells = 0

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate RMSE: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    !mean = GetMean (grid, maskReal = maskReal)

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
  
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO

ELSE
    
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

rmse = ( rmse / countCells ) ** 0.5

IF (PRESENT(nrmse)) THEN
  IF (nrmse) THEN
    mean = GetMean (grid1)
    rmse = rmse / mean
  END IF
END IF

RETURN

END FUNCTION GetRMSEfloat